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Abstract 

Based on the special properties of Liouville eigenoperators a perturbation the- 
ory for the partition sum is given. It is applicable for any temperature and 
includes the case of degenerate Hamiltonians. To demonstrate the realibility 
of the method, the second order correction to the atomic limit grand canonical 
potential of the Hubbard model is calculated and compared to results known 
from the literature. 



1 Introduction 

From the early beginning of statistical thermodynamics there was no doubt that the 
partition sum is the central quantity for describing the equilibrium quantities in physics. 
In spite of its central role there is no easy to handle perturbation theory for arbitrary 
temperature up to now. Of course, from the first days of quantum mechanics there were 
attempts in this direction. They all suffer from the non-commutativity of the perturbation 
with the unperturbed Hamiltonian, which makes it necessary to introduce an imaginary 
time via the Feynmann time ordering trick to split the exponential contained in the 
partition sum, resulting in various graph schemes for perturbation series or two time 
Green functions. Of course, perturbation theories are widely employed in all fields of 
contemporary physics and, therefore, it is hopeless to review all the developments, so 
I restrict to a brief review of the state of art in the context of the Hubbard model 0, 
^, a basic model in condensed matter physics due to its relevance for strong electron 
correlation phenomena like itinerant magnetism or, for meanwhile ten years, high T^ 
superconductivity 0. In 1980 Kubo @] published a high temperature expansion, and the 
nowadays developed cumulant expansions ||, Q and in future the incremental method 
seems to have the potential for a break through in the direction towards arbitrary 
temperatures. Nevertheless, at the moment the latter theories are mainly elaborated 
for the ground state properties, whereas the finite temperature business was left aside. 
Another way, which was shown to be equivalent to the cumulant technique ^ at least for 
the ground state, is the coupled cluster expansion [|1^. An extension to finite temperature 
was given in In p| an algebraic approach to operator perturbation theory was 

presented for zero temperature, nevertheless an extension to finite temperatures seems to 
be possible also in this line. The series expansion for the thermodynamical potential can 



2 



be generated from series expansion of the one particle green function, what of course also 
introduces (imaginary) time variables. In the following we will develope a series expansion 
avoiding this difficulty. 



2 The perturbation series for the partition sum 

We assume, that the system under consideration is described by a Hamiltonian which 
can be splitted into two parts, i.e. 

H^H^ + H^, (1) 

The partitition sum is 

Z = 5p{e-/^(^o + i?i)| (2) 
This can be rewritten in the following form 

Z = 5p|e-/^^05(^)| (3) 
with 

S{(5) = e^^Oe-^^ (4) 

The perturbation expansion for the operator S{pi) given in various textbooks on quantum 
statistics is 

= E(-l)" / d(3i d(32--- d(3nH,{(3,)H,{(32) ■ ■■H,{l3^) (5) 

n^O Jo Jo Jo 

with Hi{t) being 

It is this (imaginary) time dependence and the multiple time integrals which makes the 
calculation of the perturbation scries involved. In the following we show that one can get 
rid of this problem. For that, we consider both parts of the Hamiltonian as elements of 
the same operator space, and, therefore. Hi can be expressed in terms of eigenoperators 
of the Liouvillian belonging to Hq. The Liouvillian of a Hamiltonian is defined as 

LA = [H,A] (7) 

and an eigenoperator A fuUfiUs the eigenvalue equation 

LA = XaA. (8) 

The product of two eigenoperators A and B is also an eigenoperator, with 

L{AB) - {LA)B + A{LB)^XaAB + AXbB^{Xa + Xb){AB) (9) 



Next, we expand the S from eq (W) into a Taylor series 



oo /On 
n=0 "'■ 



Qn 



dp' 



(10) 



J/3=0 

For the first derivative of S we find 

^S{(3) = an(/5)5(/3) (11) 
with 

anif3) = e^^Oane~/^^o = -eP^^H^e~P^^ (12) 
The operator an is first order in the perturbation. The second derivative is 

= (a2i(/5) + a22(/9))5(/5) (13) 



;S{P) = — an(/3) + an(/3)an(/?) 5(/3) 



Since the derivation of an a-operator with respect to (3 does not increase the order of the 
perturbation a2i(/3) remains first order whereas 0-22 (/3) = a.ii(/5) is second order. The 
next derivative is 

S{(3) = f4«2i(/9) + 4«22(/3) + (a2i(/3) + a22(/5))an(/5)l5(/5) (14) 



d(3^ '"^^ \d(3 d(3 
Collecting together terms of the same order in the perturbation yields 

-^S{(3) = {a,,{f3) + a,,{f3) + a,,mS{(3) (15) 

Proceeding in this way we can sort the summands contributing to the expansion of S in 
a table. Since the derivatives have to be taken at infinite temperature (i.e. at /5 = 0) we 
have So = 1. Furthermore, the derivative of an a„m with respect to (3 taken at (3 = 
is nothing but applying the Liouvillian to that operator. In table 1 the column index 
gives the order in the perturbation whereas the row index gives the order in /3. Thus by 
resorting the sum one finds 

S{P) = l + S,{P) + S2iP) + SM + --- (16) 
with 

00 on 

Sm{(3) = E ^^nm (17) 

n=m 

It is easy to show that for the operators the following recursion relations hold 

O'mm = "ll (18) 

a„i = V^-'au (19) 

ail (20) 



1. 



2. 



3. 



4. 



1! 



2! 



3! 



4! 



an 
I \ 

I \ I \ 

La Xaii Lo Xaii 

C*32 Ci33 

I \ I \ r 

I XI XI 

0^41 O42 0-43, 



X 



Ct44 



Table 1: Operators a„m of the same order in the pertur- 
bation are depicted in the columns. The rows give the 
operators contributing to the same power of the temper- 
ature. The arrows symbolize the recursive connection 
between the operators. A vertical arrow stands for ap- 
plication of Lq and a slanted arrow for post-multiplying 
with ail. 



By help of the above equations (|20|) the following formula results 

n— 1— m 



a 



(21) 



fc=0 



The recursion relations allow to find closed formula for the m-th order contribution to 
the partition sum. We will show this explicit for the second order contribution and the 
third order contribution respectively. Summing the second and the third column in table 
1 respectively yields for the operators ^2 and ^3 



^2 n\ 



^3 ^' 



:22) 



Iterating eq (|2l|) one gets for the operators a„2 

( n-3 

0,n2 = Sn,2Q'22 + Sn>2 I Lq 0-22 + (a„-i_fc lan) 

I k=0 

( n— 3 



^n,2Q'22 + ^n>2 \ Lq ^022 + ^ -Lq 
I k=0 



'0 



aiijaii 



(23) 



Here the shortcut (5„>m says that n has to be greater than m. For the third order results 

O-nS = 5n,Z a33 + '^n>3 -^'0 ^ O,?,?, + 
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{n-4 / n-4-A: \ ^ 

J2 Lo ^o"'"'«22 + E 4 [(l>r'-'-'an)an] an (24) 
A:=0 \ /=0 /J 

One finds with eqs (|3|,p!7D for tlie m-tli order contribution to tlie partition sum 

ZJZ, = Sp\e-PH^Y.^^nm\=Y.^^^nm) (25) 
t n=m J n=m 

wliat can be rewritten by lielp of eq (|21|) 

)) (26) 

n=m ^- t k=0 ) 

Since for every operator A holds 

= 5rUA) (27) 

we find 

oo pn 

Zm/Zo = riO-mm) + /, — T (an-l.m-lOn ) (28) 

the next recursion step yields 

pm oo pn 

Zth/Zq = rio-mm) + /, — r((-^0 ^^m-l,m-l)o,ll) 

n=m+l ^! 
oo on n—l—m 

+ E ^ ((-LS(a„_2-fc,m-2aii))aii) (29) 

n=m+l k=0 

The recursion ends when the operators are reduced to products of an. Thus we get for 
the second order from eq ( pS] ) 



p2 oo pn 

Z2/Z0 = -Tr{a22) + E — j-(an-iian) 

^- n=3 

/32 00 <Qn 

= ^(«n) + Err((4""Vi)an) (30) 
and for the third order 

o3 00 /On 00 an n— 4 

^3/^0 = ^(a33) + Ezr((4""'«22)an) + E^E((^o(io-'-V)an)an) 

n=4 n=4 1=0 

o3 oo oo /On n— 4 

= ^(«n) + E ?r((4"'«n) «n) + E E((^o (Lr'-'au) an) an)(31) 

n=4 n=4 1=0 

Next, the perturbation Hi is decomposed into eigenoperators of Lq. For that we use the 
natural basis of the Liouville space formed by the dyades constructed from the eigenstates 
of Hq. In this representation Hq is diagonal and Hi is chosen purely offdiagonal, what is 
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always possible. In case that Hi contains components which are diagonal they are added 

to Hq. 

Hq = J2^iyW){^\ '^ith Hq\u) = e^\u) . (32) 
and 

^i = EEVI/^)(H (33) 

The basis operators |yu)(z^| are eigenoperators of the Liouvillian following 

LQ\y>i{y\ = A^^|/i)(z/| with = e^-e^ (34) 

For the operator an we get 

"11 = EE^A''^I/^)('^I ^i*^ A^i, = -Vf,u (35) 



By insertion of the above expression into eqs ( |30| , |3lD and using the linearity of the Liou- 
villian the result for the second order contribution to the partition sum is 

^2 = EEe~^^'^VK^/2(Ai) (36) 
with 

/.(AO = f + ^ (37) 
and the third order contribution reads as 

^^ V ct 

with 

^ ^ 1 ^3(Al) ^ 1 ^3(A 

3! A2 — Ai Ai — A2 A 

Here we introduced the functions i?m(A), which are determined according to 

. _ I; (^ 

n=0 

Furthermore we abbreviated 

Ai = A^^. (41) 

A2 = A^^j/ + \yyl (42) 



^3 = -E EEe"^^^' VKaW3(Al,A2) (3^ 



/3(Ai,A2) = ^ + + (39) 
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In the same way we proceeded to higher orders in the perturbation. We calculated the 
coefficients to the eighth order terms, what one finds for the m*^-order contribution to 
the partition sum is 



Zm = {-ir e-Z^^/^i V,,,,V,,,, ■ ■ ■ V,^_,,^V,^,JMu . . . , A„_i) (43) 

with 



i=l iij^A^j - ^i) \ 



and 

Afc = A^j^2 + A^2/i3 + ■ ■ ■ + ^t^k^^k+l ~ ^Mi ~ ^Mfc (^^) 

We want to emphasize, that the functions /m(Ai, ■ ■ • , Xm-i) remain finite for any number of 
vanishing A,, and, therefore, the case that the unperturbed Hamiltonian has a degenerated 
spectrum is included. Furthermore, it remains finite for all values of (3 . This becomes 
immediately clear, if one looks at the structure of the operators a„m; which do not contain 
any denominators, e.g. look at a„2 and a„3 as given in eqs (p3D and (p^) respectively. 
Eq (|43D provides a compact form for the perturbation series of the partition sum which 
can be tested easily for simple systems. As will be shown in the following the functions 
/m(Ai, ■ ■ ■ , Am-i) are nothing but the result of the m-1 integrations necessary in standard 
perturbation theory. Thus we did this integrals to infinite order. Starting from eq (^) we 
can insert the perturbation in the form given in eq (|33|) , what yields 



13 /3i /3™-i 



due to the eigenvalue equation of the unperturbed Hamiltonian for the partition sum. 
It is obvious that the result of the m integrations have to be the same as the functions 
introduced in eq (^^. Taking into account the definitions of the Aj one finds 



£ur - £u2 = Ai (47) 
'u, - £u,+i = + \ with 1 <i <m (48) 



and, therefore. 



/3lAi /32( — Ai + A2) . . . f,Pm{ — ^m-l] 



These integrals were calculated by help of symbolic computer programs. Up to order five 
(2^ terms), it was possible to bring the results symbolically to the forms given in eq. (^) 
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by help of the tools, the symbolic computer languages provide. For the sixth, seventh, 
and eighth order term we proved it by numerical calculation, for arbitrary sets of A's, 
and found reasonable coincidence within the numerical error. A different way to do the 
integrals is by help of Laplace's transformation [^. This results in a form of fm which 
seems at the first glance different. Nevertheless it can be shown to be equivalent to eq 
(H). The calculation is given in the appendix D. For comparison with the literature, where 
results are given for the series expansion of the grand potential, we have to express the 
perturbation series for the grand potential F via the perturbation series of the partition 
sum. One finds easily 

pF = (3Fo + (3F2 + (3F^ + (3Fi... 

(50) 



Z2 




(z. 


Zl 


^0 


Zo 




Zl 



This expression may be compared to the usual perturbation theory for the grand potential 



m 



f3F = /3Fo-^(-l)" Tt/A rrf/?2--- r""rf/5„(/fi(/3i)/fi(/32)---ifi(/?0)^ (51) 

„=2 -^0 

here (. . .)'^ indicates that cumulants have to be calculated. Due to the linearity of 
cumulants we can separate the r-integrations what yields together with eq 



PF = /5Fo-^(-l)'^5:---5:e-/^^/^i/4Ai,...,A„_i)x (52) 

n=2 Ml Mn 

^(^MlM2^2M3 ■ ■ ■ ^n-lMn^riMl) 

Thus it is sufficient to calculate cumulants, i.e. the linked graphs. 

3 Calculating the expectation values 

What remains to calculate is the expectation value of eigenoperator products, being of 
the form 

{AB) = ^Sp\e~f^^<^AB] (53) 

Here A is an arbitrary eigenoperator of Lq with eigenvalue and B any other operator, 
(■ ■ ■) is again the expectation value with respect to the unperturbed Hamiltonian . Since 
cyclic permutations under a trace do not alter the expectation value we find 







1 _ e-f^^^Sp < 





Sp |e~/^^o[B, A] + | - e-f^^^Sp |e-^^OAs| (54) 
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Solving for (AB) yields 

{AB) = ] {[A,BU) (55) 
If B is an eigenoperator instead of A one finds in the same way 

(AB) = -{[A,BU) (56) 

Of course one can get these results also via the standard Green function technique which 
becomes extremly simple for eigenoperators. This is shown in appendix A. If both A and 
B are eigenoperators, than one has 

Spf^e-f^^oAB^ = Sp^^e-f^^oABef^^Oe-f3Ho^^ (57) 
= Spi^ (^e-f^^oAB^ e^/^^oj = e'f^^^^ + ^b) Sp je-^-^o Ab| 

It follows that either the expectation values (AB) and (BA) vanish or the equation 

\b = — Aa for A, B eigenoperators of Lq (58) 

holds. In some models, e.g. within the Hubbard model, the operators A and B may be 
both fermionic and bosonic. If both operators are fermionic, than the anticommutator 
is suitable, since the number of operators will be reduced. In case that at least one of 
the operators A or S is bosonic, than its more convenient to work with the commutator. 
The related formula can be derived easily from eq (|55D to be 



Introducing a "parity function" P{A\B) being an odd integer if both operators are 



fermionic and even else one can unite eqs (BR 59) in the following way 



(AB) = ^(AS- (-l)^(^l^)BA) (60) 

In the above form the number of operators will be always reduced. It is of some interest 
to discuss the case when the denominator in eq (|59D vanishes. This may happen if the 
temperature goes to infinity. Since there is no reason that the expectation value (AB) 
becomes infinite, we have to conclude, that the commutator becomes zero, what represents 
the classical limit. More interesting is the case that the eigenvalue Xa is zero. But this 
implies the following statement: If the operator A is an eigenoperator to the eigenvalue 
zero, then it is diagonal or at least nondiagonal only between degenerated states. Indeed 
from expanding the operator A in the basis of Lq and applying the eigenvalue equation 
follows 

LoA = XoEEVI/^)(^l = EEV(^/.-^-)l/^)(^l = (61) 
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Taking the matrix elements yields 

= A^,{e^-e,) (62) 

what proves the statement. Of course the same holds if B belongs to the eigenvalue zero. 
Furthermore, from eq (|58| ) one can derive the following statement: If both operators A and 
B respectively are eigenoperators and at least one of them has the eigenvalue zero, then 
the expectation value vanishes or both eigenvalues are zero. Thus, for all nonvanishing 
expectation values both operators have to be diagonal, if one is. In the latter case the two 
operators commute. One can also prove the following statements, in some sense reverse 
to the above said. If two eigenoperators A and B anticommute their expectation value 
(AB) vanishes. This follows from equation (|59|) 

but this means (AB) = 0. Furthermore, if two eigenoperators A and B commute their 
expectation value (AB) either vanishes or both eigenvalues have to be zero. Indeed, we 
get from eq (|55D 



{AB) = \ {[A,BU) = ^ {AB) (64) 

Again one possible solution is {AB) = 0, but if we demand that {AB) ^ then neces- 
sarily Aa = holds. Besides Xa also has to vanish which follows from eq (|58|). Thus 
the number of operators inside an expectation value can be reduced by help of eq (|^) till 
all remaining operators are commuting and diagonal. For illustration, let us assume for 
the moment both A and B to be simple basis operators of the form 

A = \fi){fi'\ with XA = £f,-£f,' (65) 

B = \i'){i''\ with Xb = Su — Eu' (66) 

then we find 

{[ABU) = [e-P^. + e-P^'^)^-^ (67) 

(AB) = \ +; ^^ = ^W (68) 

pP\E^ bu) _|_ X ^0 ^0 

For arbitrary non-diagonal operators 

A = Y,A^,\li){p\ and B = B^.,,\iJi'){u'\ (69) 

one finds 

{AB) = Y^-^A,,B,, (70) 

If we take both A and B to be the perturbation we get together with eq (p7| ) the result 
for Z2 in accordance with eq (|36|). In appendix C we show how eq (|60D together with the 
statements given above can be utilized to evaluate expectation values within the Hubbard 
model in a systematic manner. 
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4 The one band Hubbard model 

We will demonstrate it for the Hubbard model, being the most simple lattice fermion 
model taking into account electron electron interaction. In the context of strong electron 
correlation especially the perturbation expansion around the atomic limit is of interest. 
Therefore, we focus the discussion to that case. The Hamiltonian of the Hubbard model 
for a grand canonical ensemble is 

(71) 



H = 




with 




Ho = 




Hi = 









(72) 
(73) 

Here cf^ and are the creation and destruction operators in Wannier representation. 
The chemical potential fi and the magnetic field in z-direction are introduced to take the 
effects of doping and applying external magnetic fields into account. The model has two 
exact solveable limits, i.e. the band limit with U = and the atomic limit, where t = 
holds. Since we are interested in the large U limit, we use as unperturbed Hamiltonian 
Hq, the so called atomic limit of the Hubbard model. Within this limit the electrons 
are at N independent lattice sites and the partition sum factorizes. The eigenoperators 



of the related Liouvillian Lq are the so called Hubbard operators |^ and products of 
them. Since the multisite Hilbertspace and also the related Liouville space are the direct 
product of the N single site spaces, it is enough to restrict the discussion to one lattice 
site indexed by i for the moment. For one lattice site the atomic limit Hamiltonian is 

Hi = Uui^iLiii - hi{ni^ - nil) - f^i'^i^ + "a) (74) 
and its Hilbertspace is spanned by the four eigenstates: 

\i, 0) if the lattice site is empty, (75) 

\i,u) if the lattice site is occupied with a spin-up electron, (76) 

\i, d) if the lattice site is occupied with a spin-down electron, (77) 

|i, 2) if the lattice site is occupied with two electrons. . (78) 

The related eigenvalue equations are 

Hi\i, ^) = ei^\i, fi) with /i, i/ G {0, m, c?, 2} (79) 
and 

eio = , Eia- = -ahi - /x , = -2/i + U (80) 
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From these states one can construct the natural basis of the related operator space. This 
was also first done by Hubbard |^ indicating the basis operators by = |«, /x) (z, z/|. 
The eigenvalue equation for the basis operators is 

L,xr = Arxr with ^- = e,^-e,, . (81) 

For a detailed discussion of the physics of these operators and a related diagrammatic 
technique I refer to the book of Isjumov and Skrjabin and the references therein. 
The fermion creation and destruction operators may be expressed via the basis operators 
according to 

= XO'^ + aXr-^ (82) 
c+ = Xf + aX^-'^ (83) 
n,^ = Xr + Xf (84) 

and for the atomic limit Hamiltonian playing the role of Hq in the perturbation theory, 
we get 

^0 = E((|-/^)^f -(/^ + ^^^)^r) (85) 

The related partition sum is 

Zo = l[zi with ^, = l + e/^(/^ + ^^)+e/^(/^-^^)+e/^(2/i-?7) (gg) 

i 

The partition sum factorizes into the product of N single site partition sums. By help of 
eqs (|8^j83| ) the hopping part may also be rewritten in terms of single site basis operators 

ifi = t E (Xf + aXf °Xr-2 + aX^'X^'' + X^-'^X""') (87) 
Using the eigenvalue equation for the X-operators yields 

LoXrxf"' = (Af + Af "')Xrx;'''' (88) 

Therefore, the operator an defined in eq ([T^p takes the following form 

«ii =-i EE^I.. (89) 

Here I used the abbreviations 

A^, = Xr°Xf (90) 
= aXfXr-2 ^g^^ 

Al^ = aXf-X]^ (92) 
Aj, = Xf-^Xr'^' (93) 
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What remains, is to calculate the matrix elements of an. The straight forward way is to 
define an eigenstate of Hq by help of the X-operators in the following way 

N 

m) = n^rVac) (94) 
1=1 

and to calculate the matrix elements. This lengthy calculation is interesting from a 
pedagogical point of view and we shift it to the appendix C. Here we adopt another way, 
starting from eq (PD|). Insertion of operator an in the form given in eq ( P^ ) yields for the 
Hubbard model 

r i^ja r' k^la' \ n=3 / 

= *'E E E E hiX,^{^,.K.') (95) 

r i^ja r' killer' 

The third order term becomes with eq (^Tj) 

z,iz, = -t'T^ET. E i: E MK,.,K,. + >^f.'m,.K<r'^'f'.'') (96) 

r ijtja r' i'+j'a' r" i"jtj"a" 

Proceeding in the same way as before yields for contribution of order t"^ 

ZjZo = i-trJ2---J2fMu---.>^m-i){A,,---A,J (97) 

Here Xi abbreviates the set of indices r, i, j, and a of the operator A^j^. The meaning of 
Ai, A2, and so on has a little bit changed with respect to eq (^), i.e. 

Al = Xxi = Kja 
A2 = ^Xi + Xx2 



A. = E^x. (98) 

i=l 

The remaining task is to calculate the expectation values in eq (P7|). Since all the X- 
operators appearing in the perturbation are nondiagonal, all expectation values containing 
unpaired operators vanish. Thus we have to take into account all possible systems of 
paired X-operators, where a factor -1 has to be included, if the number of commutations 
necessary to make all pairing X operators neighbours is odd, a factor +1 otherwise. This 
is nothing but Wick's theorem. Of course most of these remaining expectation values 
vanish also, since not every X-operator fits to each other and the pairing of two operators 
may be on-site non-diagonal. Nevertheless after the first pairing step every expectation 
value containing 2m, what is two times the order in t, X-operators disintegrates into a 
finite series of expectation values containing m X-operators. In case that the lattice site 
indices of all remaining X-operators are different one from each other, we have to take 
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into account all terms containing diagonal X-operators. In case that some of the lattice 
site indices are equal we have to contract them again. This way one can systematically 
find all nonvanishing contributions to an m-order expectation value. The second order 
expectation value is 



(^[,.^w) = (Xr^'Xf^X^-^-X^^""") (99) 
Due to the summation restriction i ^ j and m ^ n two contributions remain 

(xf'"^xf^^x;^™^-x^"<) = (100) 



e ^ e 



^3 



Insertion into eq (^) shows, that only four out of the sixteen terms survive. What we 
get finally is 

g Pi^^ia ~l~ ^J — it) 

H /2 i^ia — £iO + £j-a — £j2) 

Z'lZj 

H /2 i^i2 — £i-o- + £jO — £ja) 

ZiZj 

^-I3{ei2 + Sj-a) \ 

H /2 {ei2 - ei-„ + Ej^^ - ej2) (101) 

ZiZj J 

With eqs (^) we get therefore 

Z2/Z0 = t'EE f2{-cr{h,-h,)) + f^^-a{h,-h,)-U) 

„/3(2/i - U) ^f3{3fi -U- ah,) \ 
+ /2 {U - a{h, - h,)) + /2 (-a(/i, - h,)) (102) 

ZiZj ZiZj j 

Although the above formula holds for arbitrary magnetic fields, we restrict here to the 
most discussed ferromagnetic and antiferromagnetic cases and admit nearest neighbour 
hopping only. For a homogenous magnetic field, i.e. hi = hj = h, we find 

„/J(2fi - U) -,8(3fi -U- ah) \ 
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Here N is the number of lattice sites and g the number of nearest neighbours. For a 
staggered magnetic field, i.e. hi = —hj with hi = hg on the A sublattice and hi = —h^ on 
the B sublattice, we have 

g/3(2/i - U) .p{3fi -U- Xh) 
+ ^^/.(^-2AM + /.(-2AM 1 (104) 

with 

^ _ f 1 for i G A and a = +1 or for i G B and a = — 1 fl05) 
1 —1 for i G A and a = —1 or for i G B and a = +1 



We compared our second order result given in eq (|102|) to that given in |T^, |T^, |T8|, 
and found it in complete accordance. Here we restrict ourself to the second order since we 
believe it is enough to demonstrate the reliability of the presented perturbation theory. 
The calculation of higher order terms and discussion of the physics contained therein we 
shift to a forthcoming paper. 



5 Conclusion 

The main result of this paper is condensed in eq (^4]), since it contains the time integra- 
tions. We showed here the derivation via recursive relations, since this way it is palpably 
that the degenerate case is included and does not generate any problems. Furthermore, 
the recursive relation for the operators a„m can be solved easily in symbolic manner by 
help of symbolic computer algebra programs. Once knowing the factor fm, the remain- 
ing task, i.e. the calculation of the diagonal matrix elements, is straight forward. Since 
linked graphs have to be evaluated only, the further calculation steps are very likely to 
the coupled cluster expansion. In [|l^] the authors did the time integrals by symbolic 



computation. In our theory this step is economized. At the first glance one might think 
that this benefit is paid by the disadvantage that one cannot permute the indices due to 
the weigth factors /^(Ai, . . . , Am-i). This is not the case due to the inherent symmetries 
of the functions fm- Furthermore, since these functions remain finite for arbitrary sets 
of energies, the case of degeneration is included what is important for the typical models 
of strong electron correlation. Nevertheless, it remains cumbersome enough to calculate 
all the linked diagrams in higher orders. The reduction via commutations as shown in 
the appendix C is suited for symbolic computer algebra. In case of the Hubbard model 
the result for every cluster is a sum of products of Kronecker symbols and single-site 
expectation values of the X-operators, selecting a set of energy eigenvalues which specify 
the values of the A's in the functions fm- The form of the theory given in section 2 seems 
to exhibit slight differences compared to the variant given in section 4 for the Hubbard 
model. This stems from the fact, that in section 2 the perturbation was treated without 
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knowledge of its inner structure, whereas in section 4 we made use of this knowledge. If 
the special form is known, one can calculate V^i, as was demonstrated in appendix B, thus 
showing both variants to be identical. 
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Appendices 

Appendix A: Green function technique for eigenoperators 

We define the retarded Greens function with respect to the operator Hq in the usual way 

Q 

= -ze(t)([A(t),B]+) (106) 
6(t) is the step function, and 

A{t) = (J''^Hoj^^-itH^ ^ e^^-^OA (107) 
The Fourier-transformed equation of motion is than 

u{{A-B))^ = ([A,S]+)-((XoA; S))^ (108) 

Due to the eigenvalue equation holding for the eigenoperator A this can be solved for the 
wanted Greens function 

{{A-B))^ = ii^i^ (109) 
The expectation values (BA) one gets via the spectral theorem 
{BA) = lim-/ ^— -{((A;S)U,,-((A;B))_4 



-oo eA^^ + 

{[ABU) 

(,-(5\a + 1 
This is the same result as given in eq (|55[). 



(110) 



Appendix B: The straight way 

First we calculate the matrix elements of an. To this end we define an eigenstate of Hq 
by help of the X operators in the following way 

N 

iM) = n^rVac) (111) 

1=1 
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Here {/x} stands for a set {/^i, ...,/iAr}, summation over {/i} means N sums over the fii. 
The matrix elements to calculate are of the form 

m\xr^xf^'\M) (112) 

Insertion of a unit operator yields 

T.m\xrH{i}){{i}\xf^M) (113) 
{7} 

For the matrix element of a X-operator we find 

(Mi^ri{4) = 

(vaci (uxr) ( n f n f ri xA ivac) (iw) 

\«=1 / \/=i+l / \l'=N J \/=i-l / 

Now commuting the Xj operators to the center yields 

(Mi^ri{4) = 

Wi J V^i J 

Here we used the parity function introduced in eq (|60|) . The three X-operators at lattice 
site i implode according to 

X^^^X^^^Xr' = Xf\^^^6^^^^ (116) 

Since X^^ is bosonic, moving it to the upmost right (or left) will not change the sign. 
Furthermore it does not alter the vacuum, so that it may be omitted. What results is 

(wi^ri{4) = 

(_l)P(Mf..|o.,)(_i)P(.o|{.r.-)5^^^^^^,^^^^^^| f^^^oA (j^xr) Ivac) (117) 
Due to the orthogonality of the eigenf unctions we find finally 

m\xr'^\{-}) = (-i)^«^^-i°-H-i)^(-°i^^^'-")5.,.,5«v.n^... im 

Another way to simplify the same matrix element is to change X"'"* either to the left till 
it is the right neighbour of X°^' yielding 

(Mi^ri{4) = (-i)''"^^-^'"'"''^5M.«.'5.;..n^..^. (119) 



or to the right till it is the left neighbour of X^'^ resulting in 

({/^}i^ri{4) = (-i)^^°^°'»'^^^-"^5..«.5<..n^M.. (120) 

The different forms can be created also by applying the following rules holding for the 
parity function P{A\B) 

P(A|B) = P(-B|A) (121) 

P{A\BC) = P{A\CB) (122) 

P{A\BC) = P{A\B) + P{A\C) (123) 

P{A\BC) = P{A\BDDC) = P{A\BD) + P{A\DC) (124) 

Thus we get for the matrix element in eq ( |112| ) 

{7} ¥i I'^i 

^ (-l)-(-"m.}f.J(_i)-fee;iM-^)5^^^^^^,^^^^^,^^,5^,^^. n 5,,,, (125) 

with this result we find the matrix elements of the perturbation to be 

m\HMu}) = -tn^...EE(-i)^^^'^^'^-n-i)^^^|^^^^^')x 

x{Sf_L,a5ou, + O-5f,^25-auJi5^j05aUj + (y^fij-ahu,) (126) 

Here we used the fact that both and X^~^ are fermionic and we abbreviated both 
with /. The second order contribution to the partition sum reads now 

Z2 = t'j:j:j:j:j:j: e-^^^^^^uE^,} - e^^^) n 5.,. n ^m... x (127) 

ijkj (J rriy^n a' {^} {u} Ij^iJ 

x(_i)^(/IWf+i)(_l)^{/IM^^+')(_l)P(/IM^+i)(_l)P(/IM^+^) 

Next we carry out the {z/}-sum. All contributions where the indices i,j are distinct from 
m and n vanishes due to the fact that the X-operators within the perturbation are non- 
diagonal. What remains are the two contributions according to i=m, j=n and i=n, j=m. 
The first contribution also vanishes, what may be shown by doing the same steps we shall 
employ in the following for the calculation of the nonvanishing second case. Introducing 
the shortcut {/i} /ij for the set of /I's except /ij and Hj we find 

Z2 = E E E e-/^^({/^}A.) E E e"^^'- + '"^^Me,. + e., - + ..,) X 
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x(-i)^'{/IMf+i)(_i)^(/IMW+') 

By help of the rules given for the parity function eqs ( |1 2 1| , |1 22| , |1 23| , |1 24| ) it is easy to show 
that we get a plus sign. Expanding the products, we find that only four out of the sixteen 
terms survive 

Z2 = E E E e-/^^({/^}A.-) E E + '"^^Ue,. + e., - + X 

Summing over /Xj, fij, i/j, Uj and a' yields 
Z2 =t' E e-^^(WA.)^^x 

+ e-^(^^2 + + - e^-. - £,2)} (129) 

The sum over {yu} /j^ may be expressed by the unperturbed partition sum according to 

E e-^^({/^>A^) = ^ (130) 

W/ii ^^^^ 

with 2;, beeing the unperturbed single site partition sum. Thus we get exactly the same 
result of the second order contribution to the partition sum as was given in eq (|102|) . The 
reader might have got the impression, that the method is complicated due to the multitude 
of factors (— 1)^({---}|{---}) containing dummy lattice sites. This is a result of our aim 
to demonstrate here the straight forward character of the theory instead of modifying 
it to the special case of the Hubbard model. The straight calculation consists of firstly 
writing down the expression of Zm in dependence of the Ai, ■ ■ ■ , Am-i , secondly one has to 
calculate the matrix elements of the perturbation with respect to the unperturbed (many 
body) states, and thirdly one has to multiply the m*'* power of that matrix to the known 
function /m(Ai, ■ ■ ■ , A^-i)- This way we have not to evaluate any graphs or difficult Green 
functions, instead simple matrix multiplications have to be carried out, thereby specifying 
the A's, what is very comfortable as long as the dimension of the matrix is not to large, 
what is the case for instance in small cluster problems. In condensed matter systems, 
especially if one is interested in the limit N 00, other methods can be utilized, as was 
shown for the Hubbard model in section 4. 
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Appendix C: Systematic evaluation via commutations 

For an automated calculation Wick's theorem for the X-operators is not very convenient. 
A more systematic calculation starts from eq (§3). Again it should be demonstrated for 
the evaluation of the second order expectation values, having the form given in eq (|99D . 
In the first step all X-operators are fermionic and we get therefore 

{Xr'xfxi^'xf) = — -i {[xr'.xfxi^'xfu) 

= ~^ {S^'.ixf xr'xf) + {Xfxfxf)) 
e^-^i + 1 

— (5a',(xf' xr'xr') + s^.^ixfxz^'xr')) (131) 

eP\ + 1 

Here we used the fact i j. In the next step we have to calculate the product of three 
X-operators. Since in the general form used here, it is possible that e.g. Xf^ is either 
fermionic or bosonic, we have to take into account the parity function. As an example we 
show here the calculation of only one of these expectation values 

{xf xr'xf) = — ^-—^{[xf, xr'xf],) 

f-i)^(^fVr')A,, 



^ {6^.,{Xr'xf) + 6^,>{Xr'xf)) (132) 



g = ^_if(xf\xr'xr') 

Since we have i ^ j the operators Xr and Xj^ either commute or anticommute. It 
follows from the statements given in section 3 that the expectation value of two (anti- 
) commuting X-operators at different lattice sites (Xf"'x|^ ) is zero except that A^" = 

and = hold simultaneously, what means here that we have Ea = Sa' and = e^i or 
equivalently a = a' and ^ = In a condensed form this reads 

(Xr'xf ) = (Xrxf)5.„,% (133) 
For the remaining expectation value (X""X|^) one finds by a direct calculation using eq 



251) xr and X^ 



e f^^^'^ e ^^^^ 

{Xrxf) = (134) 



Appendix D: An alternate derivation of fm 

This alternate derivation was given by Walter John after critical reading the manuscript 
of this paper. 
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In the following we start from 



fmiPi -^1; • • • ; -^m-l) 



dpi dp. 



(135) 



Here we mentioned the parameter /5 explicitely. The Laplace transformed function is then 



F^(p;Ai,...,A™_i) = ldPe-Pf^fmiP;Xu...,Xm-i) 



Partial integration yields 



(136) 



{p; Ai, . . . , Xm-i] 



-pP 



dPe 



-P 



-{p-Xi)P 



-fmiP', Xi, . . . , Xm-l) 



(137) 



dP2e 



/?2(-Ai + A2) 



dPs 



The first term on the right hand side vanishes, if p is chosen adequately. Integrating by 
parts another time delivers 



Fm{p] X 



1) 



Xm-l) 
1 1 



P p - Xi 



dPe 



-{P-X2)P 



dpse 



P3{-X2 + A3) 



/33 



dPi 



(138) 

-1 

dp^ef^^^~^^-^^ 



Thus from repeated partial integrations one gets finally 



Fm {p; Xi, . . . , Xm-l) 



1 "^^ 1 



p2 fj^ p- Xi 

The original function one gets via back-transforming 

^ C+ioo 

fm{P;Xi,...,Xm-i) = ^ / dpeP^Fm{p;Xi,...,Xm-i) 



(139) 



(140) 



C—ioo 



Here C is a real constant larger than the maximum Xi. The integration path may be 
deformed to encircle the individual poles lying on the real axes. Thus we get from Cauchy's 
theorem 



fmiP', Xi, . . . , Xm~l) 



E 



-1 ppXi 



1 



+ 



1 xj nr^-^(A.-A,)^nr=-/(-A,) 



/ m-1 1 N 

^ i=i 



(141) 



The second term on the right hand side stems from the double pole at zero. Here we 
admitted only simple poles to get a concise form. There is no problem with the degenerate 
case, since in case that a pole has a higher order, one has to take higher derivatives during 
application of Cauchy's theorem. Via substitution of the exponential factor by help of eq 
(HO) one can show the derived formula to be equivalent to eq (H). 



22 

References 



[1] 

[2] 
[3] 
[4] 
[5] 
[6] 
[7] 
[8] 



J. Hubbard, Proc.R.Soc.London, Ser. A 236 (1963) 238 

M.C. Gutzwiller, Phys. Rev. Lett. 10 (1963) 159 

P. Fulde. Electron Correlations in Molecules and Solids. 

K. Kubo, Prog.Tlieor.Pliys. 64 (1980) 758 Springer Verlag, third edition, 1995 
R. Kubo, J. Phys. Soc. Jpn 17 (1962) 1100 

K. W. Becker and P. Fulde, Z. Phys.B 72 (1988) 423; see also § 
K. Kladko and P. Fulde, |cond-mat/9709044| 

A. W. Espinosa Mueller and A. R. Matamala Vasquez, |quant-ph/96060lT| and |quanF 



ph/9606013| . Version Nov. 6, 1997 



[9] T. Schork and P. Fulde, J. Chem. Phys. 97 (1992)9195 

[10] in nuclear physics: F. Coester and H. Kuemmel, Nucl. Phys. 17 (1960), in quantum 
chemistry: J. Cizek, J. Chem. Phys. 45 (1966) 4256 

[11] K. Elk, W. Gasser, Die Methode der Greenschen Funktionen in der Festkoerper- 
physik, Akademie Verlag, Berlin, 1979 

[12] G. Beni, P. Pincus, and D. Hone, Phys.Rev. B 8 (1983) 3389 

[13] M. Phschke, J.Stat. Phys 11 (1974) 159 

[14] B.H. Zhao, H.Q. Nie, K.Y. Zhang, K.A. Chao, and R. Minas, Phys.Rev. B 36 (1983) 
2231 

[15] W. Metzner, Phys.Rev. B 43 (1991) 8549 
[16] K. Kubo and M. Tada, Prog.Theor.Physics 69 (1983)1345 
[17] K. Kubo and M. Tada, Prog.Theor.Physics 71 (1984)479 
[18] K.K. Pan and Y.L. Wang, Phys.Rev. B 43 (1991)3706 
[19] K.K. Pan and Y.L. Wang, Phys.Rev. B 55 (1997)2981 
[20] W. John, private communication 

[21] A. A. Abrikosov, L.P. Gorkov, and I.Ye. Dzyaloshinsky, Quantum Field Theoretical 
Methods in Statistical Physics (Pergamon New York) 

[22] J. Hubbard, Proc.R.Soc.London, Ser. A 285 (1965) 542 



23 



[23] Ju.A. Isjumov and Ju.N. Skrjabin, Statistical Mechanics of ordered magnetic systems 
(in Russian), Nauka , Moscow 1997 



